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We propose a new quantum Monte Carlo algorithm to 
compute fermion ground-state properties. The ground state 
is projected from an initial wavefunction by a branching ran- 
dom walk in an over-complete basis space of Slater determi- 
nants. By constraining the determinants according to a trial 
wavefunction |\Pt), we remove the exponential decay of signal- 
to-noise ratio characteristic of the sign problem. The method 
is variational and is exact if |\Pt) is exact. We report results 
on the two-dimensional Hubbard model up to size 16x16, for 
various electron fillings and interaction strengths. 



PACS numbers: 02.70.-c, 71.10,+x, 71. 20. Ad, 71.45.Nt 



Quantum Monte Carlo (QMC) methods are often the 
only suitable tool for microscopic calculations of strongly- 
interacting many-body systems. They employ a stochas- 
tic means to sample the ground-state wavefunction or 
density matrix, from which ground-state or finite tem- 
perature properties can be computed. QMC methods 
have many different forms and are applied extensively to 
a wide range of problems in physics and chemistry. Ap- 
plications to systems of fermions, however, are plagued 
by the so-called sign problem [1,2], which arises from the 
combination of the Pauli principle and the use of ran- 
dom sampling. The common signature of the sign prob- 
lem is an exponentially vanishing signal-to-noise ratio as 
the system size or inverse temperature is increased. For 
years, this problem has significantly hindered simulations 
of many-fermion systems. 

In this Letter, we propose a new QMC algorithm 
to study fermion ground-state properties. We cast the 
projection of the ground state as a branching random 
walk (RW) in an over-complete space of Slater determi- 
nants. This approach combines important advantages of 
two existing methods, the Green's function Monte Carlo 
(GFMC) [3,4] and the auxiliary-field quant um Monte 
Carlo (AFQMC) [5,6] methods. A constrained path 
(CP) approximation is then imposed on the determi- 
nants, which requires that the overlap with a trial wave- 
function |\Pt) remain positive. The resulting method is 
free of any decay of the signal-to-noise ratio. It is varia- 
tional and is exact if |\Pt) is exact. The CP approxima- 
tion adopted here is based upon the positive projection 



approach of Fahy and Hamann [7] , but can also be viewed 
as a generalization of the fixed-node (FN) [8-10] approx- 
imation in GFMC. Test applications of the algorithm to 
the two-dimensional Hubbard model yield, for the first 
time, very accurate results (energy and various correla- 
tion functions) for large systems. 

Our general approach is independent of the form of 
the Hamiltonian, but we will use the familiar two- 
dimensional Hubbard model on a square (N = L x L) 
lattice for illustrative purposes: 



H = K 



cj„ c,- 



h.c. 



We use the imaginary-time propagator G = exp(— AtH) 
to project out the ground state \^o) from some initial 
state |*(°)). For small At, 



G : 



-ArK/2 e -ArV e -ArK/2 +0(At 3 



As in AFQMC, the factor involving interactions can be 
mapped into a one-body form by a Hubbard-Stratanovic 
(HS) transformation: exp(— AtV) = ^ x P(x)By(x). 
The auxiliary field x = {xi, X2, x^} introduces a fluc- 
tuating potential at each lattice site, and -P(x) is the 
probability density function for x. For simplicity we as- 
sume the discrete HS transformation [6,11], i.e., each Xi 
is either 1 or —1 and -P(x) is a constant, l/2 N . The prop- 
agator is then decomposed into a single-particle operator 
form: 



G = B K J2p{x)Bv{*)B k 



]TP(x)B(x). 



Here Bk = exp(— AtK/2) and any B implies B^ ® B±. 

For an initial state l^' ') not orthogonal to the ground 
state |* >, G n |*(°)) leads to |¥ ) at large n. In AFQMC, 
I 1 !?' )) is a single determinant and the path-integral 
(*(°)|G n |*(°)) = X; {x} f(x (n) ,x( n - 1 ),....,x( 1 )) is eval- 
uated by the Monte Carlo (MC) method, in which the 
n sets of auxiliary fields are sampled according to the 
overlap D. Our approach, instead, is to turn the propa- 
gation process into a RW in a space V. A point in D is 
\(f>) = \</>i)\</>i), where each \<f> a ) is obtained from the N a 
single-particle orbitals on the N lattice sites. Formally, 
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our procedure resembles GFMC [3,4], but the RW in the 
latter is in configuration space. 

The propagation is described by the iterative equation: 

= P(x)fl(x)|*( n >). (4) 

X 

The antisymmetric wavefunction |^( n )) at any stage n 
can be written as some linear combination of \<j>). In the 
MC process, |^( n )) is sampled by a finite ensemble of 
points {<Ajj™'} in V called random walkers. A stochastic 
realization of (4) is in principle straightforward: A x is 
generated randomly for each walker, and then its Slater 
determinant |<^jj. n ') is propagated with -B(x). Repeating 
this procedure for all walkers generates a new population 

{4> { k +1) } that samples |*( n+1 )). 

This naive sampling approach, however, is not prac- 
tical because of large statistical fluctuations, hence an 
importance sampling scheme [3,4] is required. We define 
an importance function Ot(4>) = (^t|^) as the over- 
lap of \<j>) with a trial wavefunction |\Pt)- Typically, 
|\Pt) is m general a linear combination of Slater deter- 
minants. Equation (4) can now be transformed into an 
iterative relation involving rather than where 
|#(0)) = O t (4>)\V(4>)). To ensure that the underly- 
ing propagation remains unchanged, -P(x) for a path 
\<j>') = B(x.)\<f>) must be modified: 

P(x) oc P{x)O t {4>')/O t {4>). (5) 

The probability -P(x) for choosing x now depends on 
both the initial and final positions of the path, but the 
iteration can still be carried out as a RW: The pop- 
ulation of walkers {<f>^} now represents l 1 ^™)). For 
each walker, a x is sampled from -P(x), |<^jj. n ') is prop- 
agated by -B(x), and then |<Ajj. n+1 ') is assigned a weight 

w[<f>^) = X^x-^( x )- The P rocess is repeated for all 
walkers in the current population, and the resulting pop- 
ulation represents |^ r ( n+1 )). We note that if an ex- 
act wavefunction is chosen as |^t), the normalization 
w(<f>^) = const, i.e., walkers will have no branching. 

In practice, -P(x) is difficult to sample directly. We 
circumvent this by sweeping^ through components of x 
one at a time: B^(x) = rii=i^( :c 0- The importance 
sampling procedure is implemented for each Xi, where 
it is easy to sample p(xi) oc O^iby (xi)<f>) / Ot(4>) and to 
compute the normalization p(xi). Walkers are 

stablized at suitable imaginary-time intervals by normal- 
izing and re-orthogonalizing the single-particle orbitals, 
as in AFQMC [6]. Schemes to control population sizes 
and reduce weight fluctuations are similar to those used 
in GFMC [4]. 

The determinant RW approach has distinct advan- 
tages over AFQMC. It is a true ground-state method 
(nAr — > oo) that can be easily carried out with efficient 



sampling techniques. We have frequently used nAr as 
large as 500, and At as small as 0.01. Also, at no extra 
cost, a better initial wavefunction l^' ') in the form of 
multiple determinants can be employed to reduce equi- 
libration time. Compared with GFMC, this approach 
automatically imposes antisymmetry by the use of deter- 
minants. It is plausible that the sign problem is reduced 
even without any approximation. Indeed, at half-filling 
or at U = 0, the approach is exact and completely free 
of the sign problem. In general, we expect that propa- 
gation with single-particle orbitals will be more effective 
than with isolated points in configuration space. Further- 
more, as we discuss below, the calculation of off-diagonal 
expectations is much easier than in GFMC. 

The most significant advantage of the determinant RW 
approach is perhaps the simple and practical implemen- 
tation of the CP approximation, which prevents the ex- 
ponential sign decay and provides a stable method. The 
sign problem occurs because of the symmetry between 
|\Po) and— |\Po) [7,12], which implies a symmetry between 
any pair of Slater determinants \<j>) and — \<j>). That is, 
the ground state resides in either half of V, separated by 
a nodal plane M (defined by (\Po|$) = 0) whose location 
cannot be specified a ■priori. In the limit of small At, 
imaginary-time paths in D are continuous. Therefore, if 
|¥ (o) (0)) is entirely on one side of M , a determinant can- 
not cross the nodal plane M without touching it. Because 
of the +/— symmetry, however, the total asymptotic con- 
tribution from any point on M is exactly zero in the path- 
integral. In other words, only paths confined in the half 
space affect |\Po)- I n a MC simulation, no knowledge of 
crossing is available; paths are sampled individually ac- 
cording to the absolute values of their weights. Except in 
cases where special symmetry prohibits the crossing of M 
(such as the particle-hole symmetry for (1) at half- filling), 
paths that touch M grow exponentially in number com- 
pared to confined paths. Therefore, asymptotically the 
signal-to-noise ratio vanishes. 

If the nodal plane M were known, we could simply 
constrain the RW in the correct half-space, and the MC 
procedure would become stable and yield exact results. 
Without precise knowledge of M , we seek an approxi- 
mate scheme with a trial nodal plane Mt- From the 
path-integral analysis above and the apparent analogy 
with FN, it is natural [7] to require that walkers main- 
tain a positive overlap with |\Pt)- However, because of 
the different basis spaces used in GFMC and CPMC, the 
physical context and implications of the FN and CP ap- 
proximations are different. In particular, the built-in an- 
tisymmetry and the over-completeness of D are expected 
to make CPMC behave rather differently. 

In their AFQMC calculations, Fahy and Hamann [7] 
applied the condition of positive overlap with a |\Pt)- 
The difficulty with their method is in the non-local nature 
of the propagation in AFQMC: any change of a HS field 
value affects all determinants that follow along the path. 
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In other words, updating of x at any imaginary-time n 
requires evaluation of the overlaps at all future times, 
and the simultaneous updating of all fields is required to 
find acceptable paths. From the path-integral picture, it 
is clear that the amount of computation increases rapidly 
with system size or imaginary time. 

With our RW realization of the propagation, updates 
and propagation are local in time. It is therefore straight- 
forward to place the constraint (vl^rl^) > on each 
walker at each time n. The constraint defines Mt accord- 
ing to |\Pt)- Any new walker that violates the constraint 
is given a zero weight and is thus discarded [13]. 

Once the RW has equilibrated, measurements can be 
made. For example, the ground-state energy Eq, = 
(*o|-ff|*T>/(*o|*T>, so the CP estimate of E is 

e c o = E^ n) )[(^ n) i ff i^>/(^ n) i*T>]/E^ n) )- 

71, fc 71, fc 

(6) 

To compute the expectation of an operator A that does 
not commute with H , \^o) must be used instead of |\Pt) 
in the estimator. This can be done using the princi- 
ple of forward walking [14] as in GFMC: We consider 
any walker <f> at time n and its descendents {^'} at 
n + r\. The paths <f> — > <f>' are distributed according to 
Ot {</>') /Ot {</>)• F° r each of these paths we can back- 
propagate |^t), i-e., operate the series of propagators on 
|\Pt) m reverse imaginary-time order. In contrast with 
forward walking in GFMC, where the computation of off- 
diagonal expectations can be difficult or even impossible, 
local and non-local operators are not distinguished here. 

Within the restricted half-space defined by Mt, the 
CP approach yields an eigenstate for H, i.e., 

ff|*o) = E o\^o)- Therefore, in this half-space E% = 
_(*S|H|*t)/(*§|*t> = <*S|ff|*o)/<*ol*o>- Since 
is an antisymmetric wavefunction, this variational esti- 
mate must have the upper bound property: Eq > Eq. 
We emphasize that, in contrast with lattice GFMC, the 
RW in space D becomes continuous as At — > 0, regard- 
less of the discrete nature of the original system. 

It was speculated in Ref. [7] that the quality of the 
approximation should improve as |\Pt) approaches |\Po)- 
We see that the physical meaning of M and the impli- 
cation of the CP constraint are in fact rather clear. For 
example, the CP results become exact if |\Pt) = l^o): 
Any \<j>) satisfying (\Po|$) = can be expanded as a 
sum of the excited states of the system. Imaginary-time 
evolution of such a determinant will remain orthogonal 
to |vPo)- Therefore, its asymptotic contribution is zero. 
From (6), we also see that Eq = Eq, and the statistical 
error disappears. 

To test the CPMC method, we calculated the energy 
and correlation functions for the two-dimensional Hub- 
bard model. As in AFQMC, results are extrapolated to 



the At = limit to remove the Trotter error from (2). 
For this purpose, calculations are carried out with at least 
three At values (e.g., 0.07, 0.1, and 0.15 for U = 4, and 
smaller values for larger U). The trial wavefunction is ei- 
ther a free-particle one or the Hartree-Fock solution from 
an initial state in which the j and j electrons are placed 
on the two sublattices respectively. The calculations re- 
quired CPU times ranging from a few minutes (4 x 4) to 
about 50 hours (16 x 16) on an IBM RS6000 590. 

In Table I, our calculated ground-state energies and 
the variational values from |\Pt) are shown together with 
exact results, for various 4x4 systems. We note that 
5 | 5 | corresponds to a closed shell case, while in 7 j 7 j 
the Fermi level falls in a degenerate set of free-particle 
states and the sign problem is quite severe in the AFQMC 
calculation. In fact, in this case it is rather difficult to 
obtain useful results with AFQMC, even at U = 4 [1]. 

A more stringent test of the algorithm is the cal- 
culation of correlation functions. Table II compares 
our results with available exact values. The structure 
factors S m and Sd are the Fourier transforms of the 
real-space magnetic and charge density correlations, and 

P(l) = a c i+i a c h°)/ N - Tne MC data were obtained 
with the forward walking and back-propagation tech- 
nique. The variational results are rather poor in most 
cases, implying that an extrapolation scheme [4] cannot 
be used. A poor |\Pt) a l so makes the forward walking less 
stable; the fact that our method works well is therefore 
even more significant. 

In Fig. 1, we show the computed Eq/N at U = 4 
for various lattice sizes and electron filling (n) = (iNT-p + 
N^/N. Also shown is available AFQMC data [17]. We 
see that our results, which are variational, are often in- 
distinguishable from the AFQMC data, which should be 
exact. It is, however, the fundamentally different be- 
haviors of the statistical errors that merit emphasis. In 
AFQMC the sign problem causes the variance to increase 
exponentially with N and nAr, while the CPMC method 
is stable and therefore exhibits power law scaling with N . 
For 10 x 10 systems, the CPMC error bars are 30-50 times 
smaller than those of AFQMC. For L > 10, the CPMC 
error bars are barely discernible, while AFQMC fails to 
yield meaningful results. Indeed, the only 12 x 12 data 
point by AFQMC not only has a large error bar, but also 
appears to lie above the CPMC value. 

We also computed the momentum distribution, pair 
correlation functions, and structure factors for systems 
shown in Fig. 1. The forward- walking remains well- 
behaved, and the statistical errors very small. For smaller 
systems (L = 6, 8, 10), various comparisons with existing 
AFQMC data indicated that the agreement was similar 
to that in Table II. For instance, the structures of the 
computed ^(k) vs. k curves, including peak positions, 
were in exact agreement with AFQMC results [17]. 

Very recently, ten Haaf et. al. [10] constructed a stan- 
dard GFMC FN approach for lattice fermions. It would 
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be instructive to compare it with CPMC. The advantages 
that the determinant RW has over standard GFMC, 
namely built-in antisymmetry, over-completeness of the 
basis space, and easier calculation of off-diagonal expec- 
tation values, should still hold. 

In conclusion, we have presented a variational, sta- 
ble QMC approach for fermions. Very accurate results 
were obtained even with a simple |\Pt)- Generalization 
to continuous HS fields is straightforward. We expect the 
method to be useful in a variety of applications, includ- 
ing continuum systems such as nuclei, atoms, molecules, 
etc. Algorithmic topics for further exploration include 
improvement of |\Pt) an d development of released-node 
[9] and interacting-walker [12] analogs. 
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TABLE I. Comparison of variational, exact diagonaliza- 
tion, and our calculated (CPMC) ground-state energies per 
site for 4x4 systems. Exact values for 5 f 5 J, systems are 
from Ref. [15], and those for 7 f 7 J, are from Ref. [16]. 





U 


Variational 


CPMC 


Exact 


5 T 5 i 


4 


-1.1088 


-1.2238(6) 


-1.2238 


5 T 5 i 


8 


-0.7188 


-1.0925(7) 


-1.0944 


7 t 7 I 


4 


-0.8669 


-0.9831(6) 


-0.9838 


7 t 7 I 


8 


-0.592 


-0.728(3) 


-0.742 


7 t 7 I 


12 


-0.474 


-0.606(5) 


-0.628 



TABLE II. Comparison of selected values of correlation functions for 4x4 systems. S m and Sd are the magnetic and density 
structure factors, respectively, and p is the one-body density matrix. Exact diagonalization results are from Ref. [15]. 







5 t 5 l,U = 4 






5 t 5 i,U = 8 




7 t T i, ?7 = 4 






P(2,l) 


S m (7T, 7r) 


5<i(7r, 7r) 


P(2,l) 


S m (7T, 7r) 


Sd(*, 1") 


p(2,l) 5 m (7T,7T) 


Sd(*, 1") 


Variational 


-0.125 


0.624 


0.624 


-0.125 


0.625 


0.625 


-0.103 4.4 


0.517 


CPMC 


-0.112(1) 


0.73(1) 


0.504(1) 


-0.092(3) 


0.76(1) 


0.440(3) 


-0.101(2) 2.9(2) 


0.428(2) 


Exact 


-0.112 


0.73 


0.506 


-0.097 


0.75 


0.443 


-0.101 2.2 


0.424 



-0.8 




1 1 1 1 




// [ 
// 


-0.9 




16x16^ 
12x12/ 




/ f 

/ r 
*° J // 


-1 .0 






/ 

A / 

// -^fL. 

// ~f 


/ r 
/ // 

// 

- r 
r 

// ~_ 






10x10 


// 


/6x6 








// 


y 


-1 .1 






// 

// 








8x8 


/y— G— 


— o CPMC ; 








// 

/ ♦ 


- - ♦ AFQMC ; 


-1 o 












0.5 0.6 0.7 0.8 0.9 1.0 

<n> 

FIG. 1. Energies per site vs. electron fillings from CPMC, together with available AFQMC data. The lines are to aid the 
eye. Curves for L = 8, 10, 12, 16 are shifted up by 0.1, 0.15, 0.2, 0.25 respectively. 
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